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We investigate the interplay between spin-orbit coupling and electron-electron interactions on the 
honeycomb lattice combining the cellular dynamical mean-field theory and its real space extension 
with analytical approaches. We provide a thorough analysis of the phase diagram and temperature 
effects at weak spin-orbit coupling. We systematically discuss the stability of the quantum spin Hall 
phase toward interactions and lattice anisotropy resulting in the plaquette-honeycomb model. We 
also show the evolution of the helical edge states characteristic of quantum spin Hall insulators as a 
function of Hubbard interaction and anisotropy. At very weak spin-orbit coupling and intermediate 
electron-electron interactions, we substantiate the existence of a quantum spin liquid phase. 
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I. INTRODUCTION 

Time-reversal invariant topological insulators^ - bulk 
insulators with metallic surfaces - are characterized by a 
Z2 invariantl^^ and cannot be adiabatically connected to 
trivial band insulator phases unless the single particle gap 
closes. While Z2 topological insulators (TIs) are robust 
against disordei^, rigorous and general results about the 
fate of TIs in the presence of prominent electron-electron 
interactions are limited ^ . Strongly correlated TIs as well 
as exotic time-reversal invariant Mott insulator phases 
have been predicted^ apart from more conventional 
magnetically ordered phases. 

By analogy to the quantum Hall effect two- 
dimensional TIs are also named quantum spin Hall 
(QSH) insulators. They were originally proposed to be 
realized in graphene^ and later also in HgTe/CdTe quan- 
tum wells where subsequent experiments measur- 
ing a quantized conductance established the field of TIs. 
They possess an odd number of pairs of time-reversal con- 
jugate counter-propagating edge states (the helical edge 
states)lffi2Ml. There are other promising proposals to 
stabilize the QSH effect in real materials: graphene en- 
dowed with heavy adatoms like indium and thalliuml^ 
and synthesized silicene was shown to exhibit a sta- 
ble QSH phase. Particularly interesting are monolayers 
or thin films of Iridium-based materials X2lr03 (X=Na 
or Li ) whi ch have been debated to possibly host QSH 
phase^^^^ since both spin orbit coupling (SOC) and 
electron-electron interactions are quite strong in such 
materials. All these systems have in common the under- 
lying honeycomb lattice where recently a gapped quan- 
tum spin liquid for intermediate interactions was found 
using quantum Monte Carlo (QMC) 20 22 

Very recent progress within ultracold atoms in tunable 
optical lattices establishes a second class of systems (be- 
side Iridium based materials) where topological interact- 
ing phases might be realized. This progress is two-fold: 

(i) tunable hexagonal lattices have been realizedl^ and 

(ii) different types of spin-orbit interactions are feasible 



nowl^. Additional electron-electron onsite interactions 
are considered as a standard tool in optical lattices All 
these achievements in such a rapidly evolving field indi- 
cate the demonstration of topological interating phases 
within cold atoms in the very near future. 

In this paper, we investigate the interplay between 
SOC and electron-electron interactions on the honey- 
comb lattice and combine two very paradigmatic models: 
to capture the non-trivial^band topology we consider the 
Kane-Mele (KM) modell^ without Rashba term and to 
describe interaction effects the Hubbard mo del, merging 
to the Kane-Mele-Hubbard model {KMR)^^^^^. 
Our goal is to combine the cellular dynamical mean-field 
theory (CDMFT)'^'^ and its real space extension with 
analytical approaches, present our phase diagram at half 
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FIG. 1: (color online). Phase diagram of KMH model within 
CDMFT, including the four phases: (i) topological band in- 
sulator (TBI), (ii) magnetically ordered spin density wave 
phase (SDW), (iii) non-magnetic insulator phase (SL), and 
(iv) semi-metal (SM) region which is shown (from right to 
left) for temperatures T = 0.025, 0.0125, and 0.005. Extrap- 
olating T ^ 0, the SM region shrinks to a line, see main text. 
All other phase boundaries are extrapolated to T = 0. Inset: 
Typical clusters as used within CDMFT. 







FIG. 2: (color online). Single particle gap Asp and magne- 
tization m as a function of [/ for A = 0. The various pan- 
els correspond to different cluster sizes used by the CDMFT 
method. A minimal cluster size C = 6 is required to find the 
non-magnetic insulator phase. 



filling including temperature effects, and thoroughly ad- 
dress the fate of helical edge states as a function of inter- 
actions and SOC coupling. Additionally, we introduce 
the plaquette-honeycomb model which illustrates that 
the QSH phase is also stable toward lattice anisotropy. 



II. KANE-MELE-HUBBARD MODEL 



The KMH Hamiltonian on the honeycomb lattice reads 



H-- 



U 2^ Ui^riii 



(1) 

where z, j, label the sites on the honeycomb lattice, Cia- 
is the electron annhilation operator, riia = c\^c^^^ t is 
the hopping integral {t is our reference energy scale, and 
hence we set t = 1), A the SOC, U the onsite interaction. 



and 



sgn[(di X dj). 



where d^/^ are the two vectors 



along the links from j to i {uij = ±1) . While the Hub- 
bard model respects SU(2) spin and Cq lattice symmetry, 
the SOC breaks the spin symmetry down to U(l) and the 
lattice to C3 while it leaves time reversal symmetry in- 
variant. Details of the used CDMFT method and its real 
space extension are provided in Appendix A. 

For large interactions, the phase is an easy-plane 
antiferromagnetl^. First, it is useful to start with 
a mean-field (MF) consideration of the TBI-SDW 
transition which can also be described within slave- 
rotor theory^ where the condensation of magnetic 
monopoles results in an XF-inst ability. At A = 0, 
Eq. ([1]) is spin rotationally invariant and a staggered 
magnetization pointing in any direction will be a good 
order parameter. For finite A, one can learn from the 
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FIG. 3: (color online). Temperature dependence of the phase 
diagram at SOC A = 0.02. Inset: Single-particle gap Asp and 
magnetization m vs. U is shown for A = 0.02 and T — 0.025. 



effective spin model (ie., the strong coupling limit 
of Eq. ([1]), cf. Ref. 8) that the magnetization lies in the 
XY plane. Hence we chose the order parameter accord- 
ingly, m = (5+) = (Sf) + ). From the Hubbard 
interaction we obtain n^^n^i = —S^S~ + (n^^ + nii)/2. 
A standard MF decomposition results in Hj ^ 

const, where the operators aka and bka are associated 
with the two sublattic es. Eventually one obtain s the MF 

spectrum = ± Y^^^(^)P"^^0^(^F"^^"^^^"k^ + const, 
and from there we find the MF equation. Note that g{k) 
is the nearest neighbor hopping and j{k) the second 
neighbor (spin-orbit) Haldane term^^. The transition 
line differs slightly from the analog calculation with the 
order parameter (Sf). Mainly, the form of reveals 
that the mean field approximation causes another mass 
term which does not compete with the SOC mass ^{k). 
Consequently, when passing from the TBI to the SDW 
phase we do not expect closing of the single particle gap. 
Indeed, within CDMFT no closing of the single particle 
gap is observed (while the gap has a local minimum at 
the transition). We also computed (Sf) which becomes 
finite at the TBI-SDW transition (see also App. A); note 
that this transition is of 3D XY universality at T = 0. 

For very weak SOC we identify a small phase at inter- 
mediate U which is non-magnetic, exhibits a finite spin 
gap, and is separated from the TBI phase by closing of 
the single-particle gap (see Fig. [3]). This phase is rem- 
inisce nt of the recently found quantum spin liquid (SL) 
phase'^^'^^', in particular, position and shape of this non- 
magnetic insulator phase essentially coincides with the 
SL phase found within QMC^^. Here, we shall mention 
that since the CDMFT method approximates the self- 
energy by restricting it to the chosen cluster (see inset of 
Fig.jl]), the correlation length of this SL phase cannot be 
inferred. Note that at least a six-site cluster is required 
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FIG. 4: Tc as function of A for TBI phase at [/ = 0. The blue 
hne is a hnear fit. 



to observe the SL phase, see Fig.|2] We have plotted sin- 
gle particle gap Agp and magnetization m for A = and 
various cluster sizes C. For C = 2 and C = 4 single 
particle gap and magnetization become finite simultane- 
ously while for C > 6 the onset of both quantities occurs 
for different U indicating the SL phase. The existence 
of such possible spin-gapped phases on the honeycomb 
lattice has also been addressed theoretically'^. 

In the phase diagram of Fig.[l] there is an addi- 
tional semi-metal (SM) phase where the Dirac cones of 
graphene are preserved. Using QMC as impurity solver, 
we systematically discuss temperature effects at weak 
SOC. We observe that the SL phase appears at very 
low temperatures (see Fig. 2) and that the SM domain 
prominently develops with temperature as the TBI gap 
is also very small when A ^ 0. For T ^ 0, the SM 
phase reduces to a line pointing along the A = axis 
which survives until relatively large interactions. From 
Renormalization Group point of view, this SM line at 
A = is known to be stable for weak /7^^. Hence, we 
show the SM region for T = 0.025, 0.0125, and 0.005 in 
Fig.[T] Remember that both the SL and TBI phases at 
A < 0.03 possess a small single-particle gap. In the inset 
of Fig. [3] we present the single particle gap and magne- 
tization for A = 0.02 and T = 0.025. The evolution of 
the phase diagram with temperature is shown in Fig.|3] 
for A = 0.02. 

We are using continuous time quantum Monte Carlo 
(CTQMC) as an impurity solver which enables us to ac- 
cess finite temperatures. This is clearly an advantage 
of the method since experiments are performed at finite 
temperatures but a careful interpretation of results is still 
needed. In Fig. [3) for instance, the SDW phase at finite 
T violating Mermin- Wagner theorem is clearly an arti- 
fact of the CDMFT method. In contrast, the stability 
of the TBI phase with increasing temperature is reliable 
and important for experimentalists. From an experimen- 
tal perspective, it may be interesting to know how Tc 
can be raised and for what parameter settings Tc is con- 
veniently large. Since the spin orbit gap for A = 0.02 is 
tiny, Tc < 0.05. Increasing A (and the bulk gap), how- 
ever, leads to an (approximately linear) increase of Tc up 
to Tc > 0.16 for A = 0.1 and /7 = 0; see Fig.|4] where the 
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FIG. 5: (color online). Double occupancy (ni^riii) as a func- 
tion of U is shown for A = 0.02 (red circles) and 0.2 (blue 
squares). The blue arrow (for A = 0.2) at [/ = 5.1 and the 
red arrow (for A = 0.02) at [/ = 3.9 mark the phase transition 
into the magnetically ordered phase. Note, the transition TBI 
(black arrow) to SM cannot be resolved. 



line is a linear fit. We corroborate that the critical tem- 
perature Tc follows the gap at T = 0. For finite values 
of a similar behavior of Tc is obtained. Note that Tc 
does not necessarily indicate a phase transition from TBI 
to metal; to clarify the precise nature of this boundary 
would require further work. 

Within CDMFT, we also access the double occupancy 
^occ = i^it^ii) = dF/dU where F is the free energy. 
In Fig. [5] we show that I^occ which measures the itiner- 
ant properties of electrons is an insightful quantity to 
detect the magnetically ordered phasef^^^^. The dou- 
ble occupancy is bounded as 1/4 > i^occ ^ 0- For 
/7 = 0, the ground state is a TBI with occupied lower 
and empty upper bands, |TBI) = Ilfc^It^^il^^' ^^^^^ 
{cLka.bkaV = T^iha.UkaV ^ud the matriccs are 
given by 

with 

a± = a±{k)=Af±g{k)[j{k)±e{k)]/\g{k)W (3) 

/?± = P±{k)=M± , 

N± = \g{k)\/^\g{kW + [^{k)±e{k)f . 

Here g{k) is the nearest neighbor hopping term, 7(/c) is 
the Haldane term (second neighbor spin orbit hopping), 
and £{k) = \^\g{k)\^^^^^y{k)^ is the single-particle Kane- 
Mele energy spectrum. Essentially, the matrices T^r con- 
tain the eigenvectors of the single-particle eigenstates. 
The Fermi level lies in the gap, sp = 0. Choosing a site 
i that belongs to the sublattice "a" we obtain, 

(TBI|n,^n,;|TBI) = 

kik2k3k4 
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|a_(fei)|2|a+(fc3)| = 



(4) 



where we used in the second last hne ( ^\^^^]^^-^^\^\}]^^^ ) = 
^k^kMzk^ as well as = 0' kAj^ka) = and 

i^ka^ka) = 1- ^ refers to the number of unit cells. Using 
the relations 



k 

we confirm that 
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\/3±{k)f = N/2 



= (TBI| Hitn^ |TBI) = - 



(5) 



(6) 



independent of A in agreement with CDMFT, see Fig.jsj 
Above we considered without loss of generality n^^r = 
al^a-^. Identical calculations performed with the bi^ 
operators which belong to the other sublattice yields of 
course the same result. In the opposite limit, ^ oo, 
we clearly find I^occ = as a fingerprint of Mott physics 
since we impose half filling. 




FIG. 6: (color online). Spectral function A{ky,uj) for a cylin- 
der geometry with armchair edges (L = 96, T = 0.05). Top: 
helical edge states of the QSH phase with A = 0.2 and [/ = 2.0 
or 3.0, respectively Bottom, left: for [/ = 5.2 > Uc (SDW 
phase) the helical edge states disappear while the single par- 
ticle gap remains finite. Bottom, right: A(ky,uj) inside the 
SL phase: no edge states cross the gap. The color bar corre- 
sponds to the intensity of A(ky,uj). 



III. HELICAL EDGE STATES 

To describe the edge states associated with QSH 
insulators, first we apply the concept of the heli- 
cal Luttinger liquid (HLL)E^. Hence, we linearize 
the spectrum around the Fermi points and switch 
from lattice operators to field operators tljji^{x) and 
i^Liix) which are right- and left-moving fields, re- 
spectively; we obtain for the non-interacting part 

Hq = vp J dx (^ilj\^^idxip — i^Li^dxi^i^^ ' Note that a 
standard single-particle (disorder) backscattering term 
^R-t^Li + which opens up a mass gap in the spin- 
less Luttinger liquid, is not allowed since the model is 
odd under time-reversal symmetry, = —1. Only 
two time-reversal invariant interactions are allowed: the 
forward scattering as well as the umklapp scattering 
~ i^R^i^R^i^L^i^Li which is not intrinsically present in 
the Hubbard model. Instead we shall include the for- 
ward interaction Hj = U J dx ^V^]^^ i^Rt^Li ^Li) • 
long as there is no magnetic order in the bulk, {Hq-\-Hi) 
can be solved exactly resorting to bosonization which re- 
sults in power-law decaying spin and charge correlations 
in the HLL. This result is also obtained through a spin 
wave analysis at the edges for weak interactions'^. It is 
worth mentioning that, in contrast, spin-spin correlation 
functions in the bulk decay very rapidly^. 

At the SDW transition, Hj turns into Hj ^ 
— Um J dxipl^^ipj^^ + h.c. with m = (i^li^i^L^)', apply- 



ing the bosonization procedure and introducing the Lut- 
tinger parameter K as usual, the Hamiltonian becomes: 



H 



-I 



dx- 



{d,ct>f + K {dJY 



C/m sin V Ai^cf) 



(7) 

Here, v is the renormalized plasmon velocity and a 
the lattice spacing, the field (j) contains both spin and 
charge degrees of freedom and dxO and (j) are conjugate 
variables^^. The Sine-Gordon term is a relevant per- 
turbation for repulsive interactions (since K <^ 2) and 
hence, through the pinning of the ^-field, the edge modes 
acquire a charge gap at the SDW transition. We can also 
check that the spin sector at the edges is described by an 
Ising order characterized by {Sf) ^ and {Sf) = 0. This 
shows that the SDW transition affects the charge sector 
of the HLL in a non-trivial way resulting in a metal- 
insulator transition of Kosterlitz-Thouless type. Note 
that the disordering-ordering transition in the spin sector 
also influences the charge properties of the edges. 

In addition, we performed CDMFT simulations on 
honeycomb cylinders with length x = (3a/4)L with both 
armchair and zigzag edges (L = number of sites in x- 
direction). By computing the spectral function A{ky,uj) 
we extract the edge state spectrum in the presence of 
interactions. In contrast to previous workl^ we treat 
the full microscopic Hamiltonian. We observe the follow- 
ing: (i) The plasmon velocity v associated with the edge 
modes slightly decreases for increasing U as expected 
from the HLL. (ii) The intensity of the spectral func- 
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FIG. 7: (color online). Spectral functions for a zigzag-ribbon 
at A = 0.05 and U = 2.0 (top) and at A = 0.2 and U = 2.0 
(bottom) are shown. 
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FIG. 8: (color online). Zero-temperature a-X phase diagram 
of the plaquette honeycomb model at [7 = 0. The a = 1 line 
corresponds to the KM model. Spectra for armchair-ribbons 
(L = 96) are shown at A = 0.15, a = 1.5 (top, QSH phase) 
and a = 0.48 (bottom, entrance of PI phase). Blue lines 
correspond to SM. 



IV. PLAQUETTE HONEYCOMB MODEL 



tion decreases with increasing U. (iii) The edge modes 
gap out when (Sf) becomes finite. In Fig.|6] we show 
exemplarily the edge modes for fixed A = 0.2, armchair 
boundary conditions, and = 2.0, 3.0, and 5.2. 

We further performed computations for the spectral 
function on zigzag-ribbons. Now, in contrast to the 
armchair-case, the periodic boundary conditions are im- 
posed in x-direction yielding kx as the good quantum 
number. For zigzag edges one extracts the more familiar 
spectrum^ of the KM model at small A (see Fig.[7[ top): 
the edge states connect the upper band of a gapped Dirac 
cone K with the lower band of the other gapped Dirac 
cone K' and vice versa. In addition, we have also shown 
the situation for stronger SOC A = 0.2 (see Fig.[7| bot- 
tom) where the bulk spectrum is relatively flat. Both 
spectral functions are computed for U = 2.0. The main 
findings (i) - (iii) obtained for armchair-ribbons also ap- 
ply for zigzag-ribbons. 

We also computed the spectral function on a cylinder 
for parameters U and A which belong to the spin liquid 
phase (see Fig.[6|. While the single particle gap is very 
small we did not find edge states inside the gap. While it 
is not known whether the true SL exhibits edge states, the 
non-magnetic insulator phase found within CDMFT is 
topologically trivial and absence of edge states expected. 



Next, we show that the QSH phase is not only stable 
towards relatively large interactions but also to lattice 
anisotropy resulting in the plaquette-honeycomb model, 
as sketched in the inset of Fig.jS] Hopping amplitudes 
which connect yellow plaquettes are = at while hop- 
ping amplitudes inside a cluster are t. Spin-orbit cou- 
plings A and A' = aX are analogously defined. In the 
absence of spin-orbit couplings, note that this type of 
anisotropy immediately opens a gap at the Dirac points. 
Mostly when a ^ 0, we check that the spectrum ex- 
hibits fiat bands associated with localized states on a 
given hexagon, hence resulting in the plaquette insulator 
(PI) in Fig.js] Also for large a, the system is an insula- 
tor (I) at the Fermi level. For finite SOC A, we find that 
the QSH phase is stable over a large region of the phase 
diagram (containing the KM model at a = 1). The gap 
closes at the phase transition between the PI/I and the 
QSH phase yielding a SM phase boundary. 

These SM lines are stable in the presence of finite 
electron-electron interactions U up to a critical [/c(a. A) 
what is of comparable size as Uc{l,X) (i.e., for the 
isotropic KMH model). Preliminary results for some pa- 
rameters indicate that for U > Uc{a 7^ 1,A) there is a 
non-magnetic insulator phase like in the a = 1 case. The 
detailed investigation of the interacting phase diagram of 
this plaquette honeycomb model is left for future studies. 



V. CONCLUSION 

To summarize, through a combination of CDMFT and 
complementary approaches, we have substantiated that 
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the QSH phase is robust toward prominent electron- 
electron interactions and lattice anisotropics. We have 
provided a quantitative analysis of the edge state prop- 
erties which shows that the ordering transition for the 
spins is also accompanied by a metal-insulator transition 
at the edges. For very weak SOC, we have also con- 
firmed the existence of a spin-gapped insulating phase 
(SL). Finally, we have introduced an anisotropic version 
of the KM model with the potential to host a rich phase 
diagram when including the Hubbard interaction. 
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Appendix A: CDMFT method 

In order to investigate edge states within CDMFT 
we use the real-space extension of the homogeneous 
CDMFT methocP^. We consider (nano-) ribbons {i.e., 
cylinder geometry) with translation symmetry in direc- 
tion along the edges (which we aligned for armchair- 
ribbons along the ?/-axis). 

Within CDMFT we map the original honeycomb lat- 
tice onto a A/'c-site effective cluster embedded in a self- 
consistent medium. The effective cluster model is ob- 
tained via an iterative procedure which can be started 
with an initial guess of the cluster self-energy T^{iuj). The 
effective medium represented by the dynamical mean- 
field, which is also known as Weiss function gdioj), is de- 
termined by the cluster self-energy T,{iuj) via the coarse- 
grained Dyson equation*^. Since the translation symme- 
try in the x-direction is broken, we shall use the real- 
space CDMFT coarse-grained Dyson equation for the 
inhomogeneous system: 



iuj -\- fi — t{ky) — Ti{iuj) 



(Al) 

Both gc and S are block diagonal matrices with block size 
Nc and dimension {{Nc xNx)x {Nc xN^)). I, J are block 
indices, Nc is the cluster size, and is the number of 
the clusters after performing the dimensional reduction, 
i.e., Nx corresponds to the height of the ribbon. 

As far as gc(i^) is determined, the impurity solver!^ 
can be used to compute the effective many-body im- 
purity problem to obtain the full Green's function Gc 



Eventually, by using Dyson equation of the cluster sys- 
tem T,{iuo) = g~^{iuj) — G~^{iuj), we recalculate the self- 
energy T^iiuo) to finish the iterative loop. The CDMFT 
self-consistent iterative loop should be repeatedly carried 
out, until the numerical convergence has been reached. 

Note that we use the spinor notation, ^"f" = {cj^,cj^} 
{i G cluster) and ^ = (^^)^ since a spontaneous trans- 
verse magnetization mixes the t and ^ spin-parts of the 
system. Therefore the Green's function in the CDMFT 
framework takes the following form. 



(*(n)vl/t(r,)) 



Gain 



■ (A2) 



From the a bove procedure we can easily obtain the 
quantites D^c^^^ and (5'^). The double occupancy 
^occ = i^it^ii) = OF/dU, as the first derivative of 
the free energy, is an important quantity in studies of 
strongly correlated systems. As Dqcc is directly related 
to the potential energy this is an indicator for the transi- 
tion order. In the following, we suppress the site in dex i 
of the electron operators. Using the notation of Eq. (A2) 
we find 



Docc = {c\c^clc^) = (T,c|(r3)ct(r2)c|(ri)c;(0)) 

= (T.4(r3)ct(r2))(T.c|(ri)c;(0)) 
-(T.4(r3)c;(0))(T.c|(n)ct(r2)) 

= Gn{0-)Gu{0-) - G;t(0-)Gn(0+) 



(A3) 



where < n < r2 < ra, and we used Wick's theorem. 

The transverse magnetic order parameter (Sf) reads 
(Sf) = ^(4^c.^ + h.c.) which can be expressed through 
the above defined Green's functions. 



{SUr = 0-)) = -{Gn(O-) + G;t(0-)} . (A4) 
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FIG. 9: (color online). (5") and (S") versus f/ for A = 0.1. 
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In Fig.[9]we show both (S'^) = (Sf) and (S^) = (5'f ) vs. U 
for A = 0.1. Since the SOC breaks the spin-rotational in- 
variance the magnetization is inside the XF-plane; con- 
sequently we find a finite value of (S^) at the transition 
into the SDW phase while (S^) remains zero. 

In the case of armchair-ribbons and PBC calculations 
we used a six-site cluster, for zigzag-ribbons an eight- 
site cluster, see inset of Fig.[T] In addition, we performed 



some PBC calculations with the eight-site cluster as well 
as with a 24-site cluster (see inset of Fig.[8|. We found 
good quantitative agreement between the results per- 
formed with the different clusters which underlines the 
reliability of the CDMFT method. But we would like to 
emphasize that a six-site cluster is needed to observe the 
spin liquid phase. For smaller clusters sizes (2 and 4) the 
spin liquid phase cannot be found as shown in Fig.|2] 
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